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Abstract. A method is presented for computing the 6-D phase-space 
density /(x, v) and its PDF v(f) in an N-body system. It is based on 
Delaunay tessellation, yielding v(f) with a fixed smoothing window over 
a wide / range, independent of the sampling resolution. It is found that 
in a gravitationally relaxed halo built by hierarchical clustering, v(f) is 
a robust power law, v(f) oc j-2.5±o.05^ over more than 4 decades in /, 
from its virial level to the current resolution limit. This is valid for halos 
of different sizes in the ACDM cosmology, indicating insensitivity to the 
initial-fluctuation power spectrum as long as the small-scale fluctuations 
were not completely suppressed. By mapping / in position space, we 
find that the high-/ contributions to v(f) come from the "cold" subhalos 
within the parent halo rather than the halo central region and its global 
spherical profile. The / in subhalos near the halo virial radius is more 
than 100 times higher than at the halo center, and it decreases gradually 
with decreasing radius. This indicates phase mixing due to mergers and 
tidal effects involving puffing up and heating. The phase-space structure 
provides a sensitive tool for studying the evolution of subhalos during the 
buildup of halos. One wishes to understand why the substructure adds 
up to the universal power law in v(f). It seems that the /~ 2 ' 5 behavior 
is related to the hierarchical clustering process and is not a general result 
of violent relaxation. 



1. Introduction 

Dark-matter halos are the basic entities in which luminous galaxies form and live. 
They dominate the gravitational potential and have a crucial role in determining 
the galaxy properties. While many of the systematic features of halo structure 
and kinematics have been revealed by TV-body simulations, the origin of these 
features is still not understood, despite the fact that they are governed by simple 
Newtonian gravity. 

The halo density profile p(r) is a typical example. It is found in the simula- 
tions to have a robust non-power-law shape (originally Navarro, Frenk & White 
1997, NFW; Power et al. 2003; Hayashi et al. 2004 and references therein), with 
a log slope of —3 at large radii, varying gradually toward —1 or even flatter at 
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small radii. The slope shows only a weak sensitivity to the cosmological model 
and the initial fluctuation power spectrum (e.g. Colin et al. 2003; Navarro et al. 
2004) , indicating that its origin is due to a robust relaxation process rather than 
specific initial consitions. In particular, violent relaxation (Lynden-Bell 1967) 
may be involved in shaping up the density profile, but we have no idea why this 
profile has the specific NFW shape. 

The properties of the velocity dispersion tensor is another puzzle. The 
velocity dispersion profile is slightly rising at small radii and slightly falling at 
large radii but is rather flat overall (Huss, Jain & Steinmetz 1999a; 1999b). The 
profile of the anisotropy parameter (5{r) indicates near isotropy at small radii 
that is developing gradually into more radial orbits at large radii (Colin et al. 
2000). For a sperical system in equilibrium, the u(r) and j3(r) are related to 
p(r) via the Jeans equation, but it is not at all clear why a(r) or (3(r) have these 
specific shapes. 

An interesting attempt to address the origin of the halo profile has been 
made by Taylor & navarro (2001), who measured a poor-man phase-space density 
profile by /TN( r ) = p( r ) /<?"(^) 3 ; and found that it displays an approximate power- 
law behavior, /tn cx r~ L87 , over more than two decades in r. Using the Jeans 
equation, they showed that this power law permits a whole family of density 
profiles, and that a limiting case of this family is a profile similar to NFW, 
but with an asymptotic slope of —0.75 as r — ► 0. This scale-free behavior of 
/tn(^) is intriguing, and it motivates further studies of halo structure by means 
of phase-space density. 

The simulations of the ACDM cosmology also reveal that the halos are bu- 
lit by a rougly self-similar hierarchical clustering process, where smaller build- 
ing blocks accrete and merge into bigger halos. At every snapshot, every halo 
contains a substructure of subhalos on top of a smooth halo component that 
has been tidally stripped from an earlier generation of substructure. Some of 
the important dynamical processes invloved in this hierarchical halo buildup are 
understood qualitatively, including dynamical friction, tidal stripping and merg- 
eres. However, a complete understanding of how these processes work in concert 
to produce the halo structure and kinematics is lacking. 

Attempts have been made to explain an inner density cusp using toy models 
of dynamical stripping and tidal effects during the halo buildup by mergers 
(e.g. Syer & White 1998; Dekel, Devor & Hetzroni 2003; Dekel et al. 2003). 
However, a similar halo density profile seems to be produced also in simulations 
where substructure has been artificially suppressed (Moore et al. 1999b; Alvarez, 
Shapiro h Martel 2002), indicating that the process responsible for the origin 
of this density profile might be a more robust feature of gravity. 

The issue of halo substructure has become timely both observationally and 
theoretically. Tidal streams associated with dwarf satellite galaxies are observed 
in the halos of the Milky Way and M31 and reveal their histories (Ibata et al. 
2001; this proceedings). Gravitational- lenses provide preliminary indicatons for 
the presence of substructure in halos at the level predicted by the ACDM sce- 
nario (Dalai Sz Kochaneck 2002). In contrast, the observed number density of 
dwarf galaxies seems to be significantly lower, thus posing a "missing dwarf prob- 
lem" (Klypin et al. 1999; Moore et al. 1999a). Also, the "angular-momentum 
problem" of disk galaxies (e.g. Navarro & Steinmetz 2000; Bullock et al. 2001a) 
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is probably associated with the evolution of substructure in halos (Mailer & 
Dekel 2002; Mailer, Dekel & Somerville 2002). While these problems neces- 
sarily invlove baryonic proceses, understanding the gravitational evolution of 
substructure is clearly a prerequisite for solving them. 

Aiming at the origin of halo structure, we report here on a first attempt 
by Arad, Dekel & Klypin (2004, ADK) to address directly the halo phase-space 
structure. The fundamental quantity in the dynamical evolution of gravitating 
systems is the full, 6D, phase-space density /(x, v), which intimately relates to 
the underlying Vlasov equation, and lies behind any (violent) relaxation process 
that gives rise to the virialized halo structure. Ideally, one would have liked 
to compute / free of assumptions regarding spherical symmetry, isotropy, or 
any kind of equilibrium, but computing densities in a 6D space is a non-trivial 
challenge. The state-of-the-art N-body simulations, with more than million par- 
ticles per halo, allow for the first time an attempt of this sort. ADK developed 
a succesful algorithm for measuring /(x, v), and studied its relevant properties 
and the associated systematic and random uncertainties. They then applyied 
this algorithm to simulated virialized halos in the ACDM cosmology, and ob- 
tained two surprising new results. First, the phase-space volume distribution of 
/ is a universal power-law, valid in all virialized halos that form by hierarchical 
clustering. Second, this power law is not directly related to the overall density 
profile, but is rather driven by the halo substructure. We thus learn that /(x, v) 
provides a useful tool for studying the hierarchical buildup of dark-matter halos 
and the evolution of substructure in them. 



2. More about / and v(f) 

A collisionless system is completely determined by the fine-grained phase-space 
density function /(x, v,t). The evolution of / is governed by the Vlasov equa- 
tion, 

dtf + v V x / - V x c& • V v / = , (1) 

with $(x) the gravitational potential, related self-consistently to /(x, v) by the 
Poisson's integral 

$(x) = —G [dx'dv {^4, ■ (2) 
J x X 

We therefore assume that a true understanding of the nature of self-gravitating 
collisionless systems must involve /(x, v) as a primary ingredient. 

While the 6D /(x, v) is hard to deal with, there is a simpler function, v(f), 
which is intimately related to /(x, v), yet is much simpler to handle. It is defined 

as 

v(f = fo) = Jdxdv8[f(x,v,t) - /o] , (3) 

such that v(f)df is the volume of phase-space occupied by phase-space elements 
whose density lies in the range (/, / + df). The Valsov equation ensures that 
every phase-space elements preserves its density along its path, as therefore v(f) 
is conserved. 

However, as the system evolves, phase-space patches of high / are stretched 
and spiral into regions with low /, and become thiner such that / is varying 
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over increasingly smaller scales. At some point one can no longer measure / 
but rather an average of it over some finite volume, referred to as the "coarse- 
grained" phase-space density, which we denote /. The coarse-grained v(f) is 
no longer conserved. In the course of a collapse or merger of a dark-matter 
halo, rapid global fluctuations of the gravitational potential re-distribute the 
energies of each phase-space element and lead to mixing. After a few global 
dynamical times, the potential fluctuations fade away, / stabilizes, and it can be 
viewed as the "physical" phase-space density of the system, since the microscopic 
fluctuations of / can no longer be measured or affect the gravitational potential. 

One obvious constraint on the final /, as an average of /, is that their 
maximum values must obey f max < /max- Since the initial v(f) is identical 
to v(f), if the initial v(f) vanishes for / > f max , then so does the final v(f). 
There are several additional constraints imposed by the given v(f) on the final 
v(f), specified by the mixing theorem (Tremaine, Henon & Lynden-Bell 1986), 
whose strength is rather limited because it only provides an integro-differential 
inequality constraint on v(f). 

In general, the same v(f) can describe for different systems. However, if the 
system is spherically symmetric and stationary, such that / is a function of the 
energy alone, /(x,v) = /(e), then there is a unique relation between v(f), /(e) 
and p(r). In particular, using dimensional arguments and the virial theorem, 
one can show that if p{r) oc r~°, then v(f ) oc / _/3 with 

18-6/3 

Therefore, the values in the range < a < 2, which are relevant for the inner 
regions of halos where / is high, correspond to a narrow range of (3 values, 
3 > (3 > 2.5. The value (5 = 2.5 corresponds to the singular isothermal sphere 
a = 2, while [3 = 2.8 corresponds to an a = 1 cusp, and (3 = 3 corresponds to a 
flat core, a = 0. 



3. Measuring v(f) in an iV-body System 

We wish to measure the 6D /(x, v) of a system represented by N particles of 
mass m each. Counting particles in uniform cells is impractical because even 
with 10 6 cells there would be only 10 cells along each axis. We therefore use an 
adaptive grid, where the cells vary in size and in shape to allow a proper reso- 
lution where needed. A particularly robust method of this type is based on the 
Delaunay Tessellation (Delaunay 1934), which has already been implemented in 
a cosmological context in 3D (Bernardeau & van de Weygaert 1996; Schaap & 
van de Weygaert 2000). 

A tessellation is the division of R d space into a complete covering of mutually 
disjoint convex polygons. In the Delaunay tessellation, every d + 1 points whose 
circumsphere [i.e., the (d— l)-dimensional sphere that passes through all of them] 
does not encompass any other point define a ci-dimensional polyhedron which 
makes a Delaunay Cell. ADK followed van de Weygaert (1994) in using the 
algorithm by Tanemura, Ogawa & Ogita (1983). An analysis of a halo with 10 6 
particles, which lasts about a week on a standard CPU, produces ~ 10 9 cells. A 
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typical particle is surrounded by ~ 7, 000 cells, made out of ~ 200 neighboring 
particles. 

The value of / at particle i is estimated by fi = {d + l)m/Vi, where Vi is 
the total volume of all the cells neighboring to particle i, d is the dimensionality 
(here d = 6), and the factor d+1 ensures mass conservation. In order to estimate 
/ at any point (x, v), one averages over the /j's of the seven particles that define 
the cell. The desired v (f) is obtained by first computing its cumulative version 



and then differentiating v(f) = —dV(f)/df. The sum in the above formula is of 
all volumes of Delaunay cells D u , whose average f v is greater than Jq. In what 
follows, we sometimes denote the Delaunay- measured / and v(f) by /del and 

■"del- 

When measuring the v(f) of a cosmological iV-body system, one expects 
two types of errors. First, the errors in the underlying / associated with errors in 
the numerical simulation itself, e.g., due to two-body relaxation effects, force es- 
timation, or time integration. A way to estimate these errors is by re-simulating 
the same system with different codes and with different sets of numerical pa- 
rameters. A systematic testing of this sort will be reported elsewhere (Arad, 
Dekel & Stoehr, in prep.). Meanwhile, ADK compared several different halos 
simulated with different resolutions and with different codes, and found that all 
the halos tested recover almost the same v(f). 

The other type of errors originate from the fact that we estimate a smooth 
/ from a finite set of particles using a specific adaptive technique. Here one may 
encounter both statistical and systematic errors. Some of these errors would 
decrease as the number of particles is increased, whereas other errors are an 
inherent part of the method. ADK introduced a simple statistical model based 
on Varonoi tesselation, which resembles the Delaunay technique and yet lends 
itself more easily to analytical treatment. The predictions of this model were 
then tested against numerical experiments with synthetic systems. 

It is found that /del at each particle is drawn from a probability distribution 
function of /dei//tme- In a typical realization with 10 6 particles, the width of 
this distribution is about one decade, describing the typical fluctuations of /del 
about /true- As the number of particles N — > oo, the shape of the distribution 
approaches an asymptotic limit, with a finite width of about one quarter of a 
decade, implying local fluctuations even in the infinite limit. 

It is then realized that fdei(/) can be viewed as a convolution of the true 
v(f) and a fixed window function, 



If v(f) is close to a power-law, the difference between v^ e \(f) and ^truc(/) is 
negligible over a large range of scales. This is a very useful feature of the 
method. 



V(f), 
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The relative statistical error in «dei(/) is proportional to 1/y/N, and can be 
approximated by 



with Vdei(/) = Jf°Vde\(f')df' and m = M/N. In practice, this means that 

when N > 10 6 , the statistical error is negligible for a very wide range of /. 
Moreover, in regions where there are large statistical errors, these errors are 
usually overwhelmed by systematic errors. 

4. A Universal Scale-Free v(f) 

The v(f) is measured from several different halos, in three different mass ranges, 
simulated within the ACDM cosmology with two different iV-body codes. We 
focus on virialized halos out to slightly outside the virial radius. The value of 
/ that can be crudely associated with the virial radius of a given halo, where 
one may expect a qualitative change in the behaviour of v(f), is estimated by 
/vir = vr -3 / 2 /9/(T 3 , with p and <t the mean virial quantities. On the other side, the 
reliable- measurement range is marked by / = /2o%> below which the statistical 
error in v(f) due to the DTFE is below 20% according to eq. 7. This formula 
was verified by the v(f) of mock systems with I0 5 particles. The statistical 
error is expected to be practically negligible in the range / v ; r < f < /2o%- When 
tested using synthetic datasets for which the true v(f) is known, we learn that 
with 10 5 particles the method recovers the true v(f) very well over a range of 
5-6 decades in /, while with 10 6 particles the range spans ~ 10 decades. 

The results described in ADK are based on three different cosmological 
simulations, two using the ART code (Kravtsov, Klypin & Khokhlov 1997), and 
the third using the TPM code (Bode, Ostriker & Xu 2000; Bode & Ostriker 
2003). The assumed cosmological model is the standard ACDM with Q m = 0.3, 
Ua = 0.7 and h = 0.7 today. Halos were sampled from the simulatied periodic 
boxes of sides L = 1, 25 and 320/i~ 1 Mpc, which we denote LI, L25 and L320. 
The simulations are by Colin et al. (2003), by Klypin et al. (2001), and by 
Wambsganss, Bode & Ostriker (2004) and Weller, Bode k Ostriker (2004). The 
force resolution is 87, 140 and 4700pc. The particle mass is 7 x 10 3 , 1.2 x 10 6 
and 2.6 x 1O 9 M0. The halo masses analyzed from these simulations correspond 
to dwarf galaxies (10 9 — 10 10 Mq), normal galaxies (~ 1O 12 M0) and clusters 
of galaxies (~ 10 15 M Q ) respectively. The ART halos were sampled by ~ 10 6 
particles within ~ l.li? v ; r , while the TPM halos were sampled by only 4 x 10 5 
particles within the corresponding radius. The LI halos were analyzed at z = 
2.33. More details about the simulations and the halo properties are summarized 
in Table 1 of ADK. The v(f) curves for these nine halos are shown in Fig. 1 

The figure shows that in every halo v(f) is well fit by a power-law, 




(7) 



v(f) oc / 



-2.50±0.05 



(8) 



over 3 to 5 decades in /. It is typically valid between about f v - u and slightly 
below /2o%- Outside this range, v(f) gradually deviates downward. In the low-/ 
regime the deviation is associated with departure from the virial regime, while 
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Figure 1. The volume distribution of phase-space density, v(f), for 
each of the nine halos analysed in ADK. The curves were shifted to 
coincde at / = /_2, where the local log slope of v(f) is —2, and were 
then shifted vertically by 4 decades relative to each other. A power-law 
line v(f) oc /~ 2 ' 5 is shown on top of each curve. Marked on each curve 
are the virial-radius level f viV and the 20% statistical error limit /2o%- 



the high-/ deviation is consistent with being due to the limited mass resolution, 
as indicated by / 2 o% and by the error analysis of ADK. The high-/ deviation 
from the power-law tends to occur at a smaller / value in L25, and even smaller 
in L320, due to the fact that M v ; r is smaller respectively. 

The power law in v(f) does not show a significant dependence on halo 
mass. There may be a marginal trend for slight steepening of v (/) as a function 
of mass, but only from steeper than /~ 2 ' 45 at ~ 10 9 Mq to flatter than /~ 2 ' 55 at 
~ 10 15 M & . This indicates relative insensitivity to the exact slope of the initial 
fluctuation power spectrum, which varies across the range from dwarf galaxies 
to clusters of galaxies. Additionally, the fact that we obtained essentially the 
same v(f) from simulations using two different numerical codes, indicates that 
the shape of v(f) is not an artifact of a particular simulation technique. 



5. Substructure 

When /(x, v) is function of the energy alone, and the halo is spherical and 
isotropic, the power-law v(f) oc /~ 2 ' 5 implies via eq. 4 that the real-space density 
profile must also be a power law, in fact an isothermal sphere p[r) oc r~ 2 , at 
least over some finite range in r. This is clearly not the case (§1), indicating that 
/ is not a function of energy alone, and the system must deviate from spherical 
symmetry or isotropy. This could be due to the clumpy substructure of the halo, 
where the surviving subhalos contribute high phase-space densities to v(f), thus 
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Figure 2. Density maps of dwarf halo Lis in a slice of thick- 
ness 0.4i? v j r . Left: real-space density. Right: phase-space den- 
sity. The units in the colour key are \og{p/[M & Mpc~' i \) and 
log(//[M0Mpc~ 3 km -3 s 3 ]) respectively. The very-high / values are 
found inside clumps which are typically far away from the halo center. 

making it shallower than expected from a smooth system with an inner density 
slope flatter than —2. 

Figure 2 shows density and phase-space density maps of an equatorial slice 
from one of the halos, in which the color represents the densities p and / averaged 
over the values assigned to each particle in a small real-space volume (of side ~ 
0.005i? v ; r ). The p of each particle was calculated using a 3D Voronoi tessellation 
(van de Weygaert 1994), which is similar in its adaptive nature to the Delaunay 
tessolation used to estimate /. While the real-space density maps are dominated 
by the familiar smooth trend of density decreasing outward, perturbed by several 
tight clumps throughout the halo, the global trend with radius becomes much 
less prominent in the phase-space density maps, with the subhaloes contributing 
the highest peaks, especially in the outer regions of the halo. Moderate / peaks 
are found everywhere (reddish yellow), and the very high peaks (bright yellow) 
are preferentially found in the periphery. The central peak in / is quite modest 
in comparison; the elongated structures near the center of the shalo are most 
likely merging subhalos. 

Figure 3 shows p and / associated with a random subset of the N-body 
particles as a function of their distance r from the halo center. A large portion 
of these particles follow the global trend of decreasing density with radius - 
they could be associated with a smooth-background component, for which / is 
approximately a function of energy alone. At radii r > 1 kpc, the high-/ values 
come in as "spikes" corresponding to the subhalos. While the spikes in p reach 
values comparable to the central peak, the spikes in / could be more than 100 
times higher, indicating that the subhalos are both compact and cold. 

The spikes seem to be lower and broader as they get closer to the halo center, 
and they completely blend into the smooth background inside ~ 2 kpc. This 
indicates that the subhalos phase-mix and lose their high phase-space densities 
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Figure 3. Densities as a function of radius for dwarf halo Lis, using 
a random set of 4% of the particles. Left: real-space density. Right: 
phase-space density. The background particles define a general trend of 
decreasing density with radius, while the spikes correspond to subhalos. 
The phase-space density spikes are higher than the central peak because 
the subhalos are cold. They become shorter and broader at smaller 
radii, indicating heating by tidal effects and mergers. 



as they approach the halo center. This seems to be the natural result of mergers 
and tidal effects, which puff up the subhaloes and especially heat them up. 

As a first attempt at trying to understand how the power law v(f) oc f~ 2 ' 5 
in a halo of mass M may originate from its substructure, one may simply add up 
the typical contributions from the general population of halos of different masses 
m smaller than M, as predicted in the ACDM cosmology. Based on cosmological 
N-body simulations (e.g., Moore et al. 1999a; Ghinga et al. 2000; De Lucia et 
al. 2004), and in accordance with the Press-Schechter (1974) approximation, the 
mass function of small-mass halos can be approximated by dn/dm oc m~ 7 , with 
7 ~ 1.8 — 2.0. Supported by the simulations, one assumes that the average 
density profiles of halos of different masses are simply scaled versions of each 
other (e.g., Navarro et al. 1996) If all the halos form at the same time, they 
have the same characteristic real-space density p m = p, so their typical radii 
scale like r m oc m 1 / 3 . Based on the virial theorem, the velocity dispersions then 
scale like a m oc m 1 / 3 . Therefore, the typical phase-space volume of a halo of 
mass m scales like V m , and its typical phase-space density is 

f m oc m/V m oc m^ 1 . With this scaling one can show that v(f) should be flatter 
than v(f) oc /-( 4 -t). For 7 > 1.8, this means that v(f) is shallower than / 2,2 , 
which is significantly shallower than the measured J- 2 - 5 ± - 05 . Including the fact 
that small halos form first in a ACDM scenario makes a negligible difference. 
One can conclude that the population of subhalos within a bigger host halo 
must be different than the general population of halos. As a result of mergers 
and tidal effects, the subhalos are expected to have a different mass function, 
their shape properties are likely to vary differently with mass, and both effects 
probably vary with radius. Thus, the phase-space density is likely to provide 
a useful tool for studying the dynamical evolution of subhalos in their parent 
halos. 
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6. Discussion 



It is important to verify that these results are not numerical artifacts. Based on 
the error analysis and tests with mock datasets, we believe that the v(f) mea- 
sured by the DTFE algorithm genuinely reflects the true phase-space properties 
of the given iV-body system over a broad range of /. The question is whether the 
phase mixing suffered by the subclumps is an artifact of numerical effects such 
as few-body relaxation, leading to underestimated inner densities and/or over- 
estimated internal velocities (Binney 2003). The apparent agreement between 
simulations run with different codes and different resolutions is encouraging. In 
order to specifically address the effect of two-body relaxation, we intend to run 
twice a simulation of the same halo with the same number of particles but with 
a different force resolution (ongoing work with F. Stoehr). 

Assuming that the simulations genuinely reflect the true physical behaviour, 
the origin of the robust power-law shape of v(f) from the merging substructure 
becomes a very interesting theoretical issue. As demonstrated in §5, a simple 
model using the mass function and the scaled profiles of the general halo popu- 
lation in the ACDM scenario does not reproduce the correct power law. This, 
and the apparent trend of the / spikes with radius, indicate that the structural 
and kinematical evolution of the subhaloes in the parent halo are important. 
Studies of tidal heating and stripping may be found useful in this modelling. 

It would be interesting to follow the phase-space evolution and the con- 
tribution to the overall v(f) by a single, highly resolved subhalo, or many of 
those, as they orbit within the parent halo and approach its center. This may 
help us understand the nature of the interaction between the parent halo and its 
subhaloes, and the origin of the v(f) power law (ongoing works with E. Hayashi 
and with B. Moore). 

We saw that the power-law behavior of v(f) is limited to the virial regime. 
It would be interesting to learn how this shape evolves in time as the halo 
virializes. A preliminary study (to be concluded and reported in another paper) 
indicates that in the intermediate-/ regime the v(f) of a pre-virialized system 
is significantly flatter than /~ 2 ' 5 , while in the high-/ regime it drops in a much 
steeper way. The /~ 2 ' 5 behavior seems to be a feature unique to virialized 
systems. 

We learned that in the haloes that are built by hierarchical clustering, the 
power-law behavior v(f) oc /~ 2 ' 5 reflects the halo substructure. It would be 
interesting to find out whether this power-law behavior actually requires sub- 
structure, or it is a more general phenomenon of virialized gravitating systems, 
valid independently of substructure. One way to answer this question would be 
to analyse simulated haloes in which all fluctuations of wavelengths smaller than 
the halo scale were removed, resulting in a smooth halo formed by monolithic 
collapse, with no apparent substructure in the final configuration. As described 
in §1, such haloes are known to still have NFW-like density profiles in real space, 
and one wonders whether they also have the magic power-law v{f). There are 
preliminary indications for a steeper v(f) in this case (Arad, Dekel & Moore, 
in preparation). If confirmed, it would indicate that the /~ 2 ' 5 behavior, while 
insensitive to the exact slope of the initial power spectrum, is unique to the 
hierarchical clustering process, and is not a general result of violent relaxation. 
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Our current results are just first hints from what seems to be a promising 
rich new tool for analysing the dynamics and structure of virialized gravitating 
systems. The analysis could become even more interesting when applied to 
haloes including the associated gaseous and stellar components. 
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